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In this paper we combine partial-wave ('modal') methods with a wavelet analysis of the CMB 
bispectrum. Our implementation exploits the advantages of both approaches to produce robust, 
reliable and efficient estimators which can constrain the amplitude of arbitrary primordial bispectra. 
This will be particularly important for upcoming surveys such as Planck. A key advantage is the 
computational efficiency of calculating the inverse covariance matrix in wavelet space, producing an 
error bar which is close to optimal. We verify the efficacy and robustness of the method by applying 
it to WMAP7 data, finding f^l = 38.4 ± 23.6 and /^l = -119.2 ± 123.6. 



I. INTRODUCTION 

In inflationary scenarios, any measurable deviation of the primordial density fluctuation from Gaussianity offers 
a window onto the underlying physics. Observable deviations require violation [Ij of at least one assumption of the 
simplest inflationary model: (i) a single, canonically normalized light degree of freedom; (ii) slow-roll dynamics, and 
(iii) the Bunch-Davies vacuum state on deep subhorizon scalesj^ For comprehensive reviews see, eg., [SHS]. 

The prospect of recovering deflnite information about inflationary microphysics has made the study of non- 
Gaussianity an area of active research. Much of this activity has been stimulated by maps produced by all-sky 
CMB experiments, including the Wilkinson Microwave Anisotropy Probe (WMAP) |6j. To maximize the value of 
these datasets requires methods capable of extracting the primordial non-Gaussian signal. The amplitude of this 
signal is conventionally expressed as /nl- 

An optimal bispectrum estimator for /nl was developed in Ref. [7] and implemented by Smith et al. For 
primordial non-Gaussianity in the local mode (to be defined in Section [TT| below) it gives the constraint /nl = 38± 21 
from 5-year WMAP data. Unfortunately, optimality of the method requires calculation and inversion of a pixel-by- 
pixel covariance matrix. Since WMAP maps contain 0(10®) pixels this is an onerous task. The calculation can be 
simplified by approximating the covariance matrix as diagonal, but the effects of anisotropic noise and masking make 
this approximation degrade with increasing resolution. The diagonal approximation is unlikely to be satisfactory for 
Planck. 

Alternative approaches exist, which aim to reduce the calculational burden at the expense of a sub-optimal error bar. 
Fergusson, Liguori & Shellard suggested that the bispectrum could be decomposed into partial waves or 'modes' [S]. 
In this approach, an efficient inverse covariance weighting has been proposed which renders the estimator closer to 
optimal [10 . But a potential drawback is the necessity to remove cross-terms (mainly generated by anisotropic noise) 
by subtracting a linear correction. For an experiment such as Planck the required cancellation may be very precise if 
the diagonal approximation is used for the inverse covariance matrix. A preliminary application of the method has 
reduced the error bar from A/nl = 29.5 to A/nl = 27.6 [TT| . 

Instead of partial waves one can use wavelets and needlets, which have a long history of application to CMB 
analysis [T^HH]- Like a decomposition into partial waves, the advantage of these methods is compression of the 
WMAP data from 0(10®) pixels into 0(10'^) wavelet coefficients, for which calculation of the inverse covariance 
matrix is comparatively trivial. Despite this, the error bars achieved by this method are only marginally worse 
than those produced by the pixel-by-pixel approach. Further, Donzelli et al. |17j have shown that, for wavelets, 
the linear correction term described above is not necessary because scale-by-scale subtraction of the mean for each 
coefficient produces effective decorrelation (see also Ref. [E]). Even in the absence of these desirable properties, 
wavelet approaches would be interesting because the use of complementary methods helps establish sensitivities to 
contaminants such as foreground noise and masking. For example, they have proven useful in performing an exacting 
noise analysis [19]. 

These properties motivate wavelet approaches to the CMB bispectrum. However, we usually wish to use the 
CMB to constrain models of early-universe inflation. These do not generate predictions for the CMB directly, and 
therefore it is not simple to compare them to wavelet coefficients recovered from the CMB: they predict the correlation 
functions of the curvature ■perturbation, C(k), which can be translated into the gravitational potential $(k). The most 
observationally important of these are the two- and three-point functions, ($(ki)$(k2)) and ($(ki)$(k2)$(k3)), and 



^ Some non-Gaussianity is always produced by post-inflationary gravitational reprocessing of the density fluctuation. In Einstein gravity 
this is expected to produce a signal of order /nl ~ C(l)i but may be larger in modified theories of gravity [2j. 
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it is the parameters of these correlation functions which we wish to estimate from the data. But to do so they must 
be converted into predictions for statistical properties of the CMB anisotropies. 

Currently, the only practical way to carry out this conversion is to systematically approximate each correlation 
function using the partial- wave expansion suggested by Fergusson et al. [5], described above. As we will explain in 



Section II B this has the effect of rendering the calculation numerically tractable. Funakoshi & Renaux-Petel ^20^ have 
recently detailed a formalism to compute a suitable expansion directly from the rules of the Schwinger (or 'in-in') 
formulation of quantum field theory. Once this has been accomplished it is straightforward to write an estimator 
for each parameter appearing in (<I'(ki)<I>(k2)<I'(k3)). We will review some aspects of the modal decomposition in 
Sectio n [lH and explain how suitable estimators can be constructed. For further details we refer to the literature [3 

[II1E1H21]. 

In this paper we take the logical step of combining a modal decomposition of each primordial n-point function with 
the use of wavelet-based CMB estimators. The CMB analysis is performed in wavelet-space and takes advantage of an 
inverse covariance matrix which is numerically much less challenging than the pixel-by-pixel case. A change-of-basis 
matrix allows us to map from wavelet-space to modal-space. This approach exploits the computational benefits of a 
wavelet-based estimator but simultaneously allows comparison to the predictions of primordial inflationary models. 

Summary. In Section |ll] we explain the partial-wave or 'modal' expansion technique and review its use in CMB 
bispectrum analysis, and in Section [III| we summarize the methodology of wavelet-based estimators. In Section [TV] we 
describe a prescription for projecting directly from the modal expansion of an arbitrary primordial bispectrum to the 
CMB bispectrum. We also review the use of the modal expansion to create simulated maps. 

In Section |V] we begin the application of partial- wave expansions to the wavelet-based estimators of Section |III[ In 
Section [VT| we describe an implementation of this prescription for the 7- year WMAP data up to /max = 1000. Con- 
straints on the amplitude of the constant, local, equilateral, flattened and orthogonal models are given in Section [VII[ 
We conclude in Section IVlIII 

Sections [TT[jIV| are a summary of the literature, and have been included to flx notation and make our presentation 
self-contained. Readers familiar with the 'modal' methodology and wavelet-based CMB analysis may wish to proceed 
directly to Section [Vj making use of references to preceding sections where necessary. 



II. CMB BISPECTRUM AND PARTIAL- WAVE TECHNIQUES 



A. CMB bispectrum 

In linear perturbation theory, the spherical harmonic transform of the CMB temperature map AT(n)/r may be 
expressed in terms of the primordial gravitational potential (f>, 

—MkMmi.ik), (1) 

where the unit vector n determines an orientation on the sky, and AT{h)/T ~ ^^^^ a/m^;m(ft)- In what follows we 
work up to /,„ax = 1000. The Ai(k) are transfer functions, which map from primordial times to the surface of last 
scattering. They are computed by solving the coUisional Boltzmann equations using publicly-available codes such as 
CAMB 125^ and CLASS [26 . 

Angular bispectrum. The CMB bispectrum is defined as the three-point correlation function of the aim, namely 
■Smim^ama = (oiimi ajamaa/gma ) ■ K Can be written 

B'A[t.^, = (47r)3(-i)'^+'=+'3 1 (^p^ ^Ai^(^k,)Y,*^^{k,)^ ($(ki)<I>(k2)$(k3)) . (2) 

We recall that the primordial two- and three-point functions satisfy 

($(ki)$(k2)) = (27r)3<5(ki + k2)P<i,(fci) (3) 
($(ki)$(k2)$(k3)) = (27r)3^(ki + k2 + k3)B$(A:i, fc2, h) . (4) 

Each inflationary model predicts a specific form for P$ and -6$. A typical model will predict the appearance of a 
finite number of momentum-dependent combinations (or 'shapes') in B^, with amplitudes that depend on parameters 
of the model. These shapes can be regarded as similar to the different Mandelstam channels in 2 — 2' scattering, 
or the structure functions of the hadronic tensor W^'' in QCD. By constructing an estimator for the amplitude of 
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each shape we obtain observational constraints on these parameters, and in some cases it may even be possible to 
rule out a model entirely. However, before any comparison with observation we must first translate Eqs. ([3|-(|4]) into 
statistical properties of the CMB anisotropy. 

The (5-functions in ([Sjl-Q enforce momentum conservation. For the bispectrum this requires that the momenta 
form a closed triangle, and implies that may be expressed as a function of the ki alone. To express the J-function 
in multipole space, we use the identity 



li'tni ^1=1 ^ i=l 



(A), (5) 



where jn{x) is a spherical Bessel function and df2 is an element of area on the sphere. After substitution into ^ we 
conclude 



(6) 



5m/4r„3 = 0) j dxdA:idA:2dfc3 (xA^ifcafcs)' (^i)A;, (fc2)A,3 (fc3)S<i,(fci, fca, fcg) 

To simphfy we note that the Gaunt integral is defined by 

Gl^[t%m, ^ J dnii) Y,:^^{±)Y,l^^{±)Y,l^^ix) = h,,i,,, ) , (7) 



where ( ] denotes the Wigner 3-j symbol, and hi,i,u ^ / (^^i + + l)(2/3 + 1) / ?i h h\ 

\ mi 1712 J V 47r \^ y 

Gaunt integral is the analog of the Dirac 5-function in multipole space, and imposes constraints on the h. Finally, we 
define the reduced bispectrum, bi^i^i^, to satisfy 

Tjhhh _ nhhh h IQ\ 

and it follows that 

biMi^i^^ y dfcidfc2dA:3 (fciA:2A;3)^Aii(fci)Az2(fc2)Az3(fc3)B$(A:i,A;2,fc3) j <lxx^ji^(kix)ji^{k2x)ji^{k2,x). (9) 

Shape function. We define the 'local' bispectrum to satisfy 

Si°^(fci,fc2,A:3) = 2(p$(fci)P$(fc2) + P*(fci)P*(fc3) + i'$(A:2)P$(A:3)) • (10) 
For any bispectrum we can define a dimensionless 'shape' function by the rule 



I. u \ - B^{ki,k2,k3) 



This choice is arbitrary: we could equally well have defined a shape function by comparison to a fiducial bispectrum 
other than B^'^. In the modal decomposition literature, the choice (A;ifc2fc3)^i3$(fci, ^2, ^3) is often made. In this 



paper we adopt (11) for numerical purposes, because it often proves more stable. To clearly distinguish our choice 



when comparing with the literature we also define a canonical shape function 5*$, 

fc2, fc3) = ^^^y fc2, fc3) , (12) 

where the normalization constant TV is adjusted to ensure S${k, k, k) = 1. With these choices the reduced bispectrum 
may be written 

Kl.l, - 6 f-) / dfcidfc2dfc3 (fclfc2A:3)'A;,(fci)A;,(fc2)A,3(fc3)P*(fcl)P*(fc2)4'"^(fcl,fc2,fc3) 

J (13) 
X / dxx'^ji^{kix)ji^{k2x)ji^{k3x). 
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B. Primordial decomposition 



Eq. ( 13 ) shows that conversion of the primordial two- and three-point functions into predictions for the sta- 
tistical properties of the CMB is simphfied whenever the shape function is separable, ie., of the form 3^^°'^^ = 
X {ki)Y {k2) Z (k^) + perms. In such cases the fci, ^2 and integrals in Eq. (13) can be decoupled, greatly reducing 
the computational time. 

To take advantage of this simplification, Fergusson, Liguori & Shellard suggested that an arbitrary (not necessarily 

separable) shape function S'^'°'^'' could be decomposed into a basis of separable partial waves [9 . The precise choice 
of basis functions is arbitrary, but should be chosen to achieve good convergence with a small number of terms. 
Fergusson et al. used a set of orthogonal polynomials qn{k) to write 

S{ki,k2,k3) =^a^rsQ{p{ki)qr{k2)qs){k3), (14) 

prs 

where bracketed indices are symmetrized with weight unity. For details of the construction of the qn{k) we refer to 
Ref. 

The physical region where the ki form a triangle corresponds to a domain V defined by 

2max(A;i, ^2, fcs) < fct , (15) 

where kt = ki + k2 + k^ is the perimeter of the momentum triangle. To simplify formulae it is convenient to introduce 
a multi-index n which runs over unique triplets (p, r, s). (By 'unique', we mean triplets which generate a unique com- 
bination q^pqrqs) after symmetrization.) Defining = q[p{ki)qr{k2)qs){k3), we write S'£°'^'' = J27i'^nQn{ki, k2, k^). 
Finally, we introduce an inner product on the physical region by the rule 

((/i5»= / dv f{ki,k2,k3)g{ki,k2,k3)uj{ki,k2,k3), (16) 



where dv is an element of area on V and w is a weight function which can be adjusted to suit our convenience. 
Using this inner product we define a matrix 7 such that 

7nm = {Qn.Qm}- (17) 

The Qn are not themselves orthogonal. Therefore, although is symmetric, it will not typically enjoy other special 
properties. But if the Qn have been chosen appropriately it will be positive-definite and invertible, in which case the 
coefficients of the separable expansion ( 14 ) can be computed, 

"? = E«4'°'^Q™))^™n- (18) 
m 

The accuracy of this expansion is limited by the number iVmax of modes used. For a shape function S and an 
approximation using N modes, a measure of convergence can be obtained by evaluating the ratio 

C(S, Sj,) ^ , ^^^'^^^^ : (19) 

For those models which have been studied in the literature, only 0(30) modes are required to achieve an accuracy of 
at least 90- 95% [9|. 



C. CMB analysis 

Given the expansion 5'^'°'^'' = J2n o^nQn, the reduced CMB bispectrum (13) becomes 

biM,=6hL J dxx^^^/'{x)qi''^'^{x)q[-'^''{x), (20) 

n— (p,r,s) 



where the summation over n is restricted to unique triplets, bracketed indices are again symmetrized with weight 
fur 

2 



unity, and the functions (jp^'*' and (jp are defined by 



g(2)'(x) = - dk k\p(- )Ai{k)ji{kx) , and q^f'^'ix) = - dk k'P„{k)qp(- )Ai{k)ji{kx) . (21) 
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The separability of the expansion reduces the integral for h^^i^i^ from four to two dimensions. 

With current experiments, the signal-to-noise available for each multipole is too weak to measure the components 
Instead, it is conventional to construct an estimator for the amplitude of each momentum 'shape' 



of bi^i^i^ directly. 



predicted by the underlying inflationary model. Any such estimator sums over many components of h-^i^i^ and can 
achieve an acceptable signal-to-noise. The optimal estimator for the amplitude of a fixed bispectrum shape bi-^i^i^ is 
proportional to [7] 



S""' = E bi.i^hGl^^.m, \(c-'a°^^) (c-'a°^^) (c-'a°^^) - 3CrX,i.^.< 



bs 



(22) 



where 'obs' indicates values recovered from observation, and C, 



-1 

hmi ,l2ni2 



denotes the inverse covariance matrix. As we 



have explained, it will typically be non-diagonal because of mode-mode coupling induced by the mask and anisotropic 
noise. If we impose the diagonal approximation discussed in the Introduction (Section [l]), the estimator £ reduces to 



Dbs ^obs „obs 



2^3 

ra2m3 



^obs o/~i n^*^^ 



(23) 



Its expectation value is (£-approx^ ^ EuiKhhKhh)ViChCi,Ci,). 

In the next section we explain the construction of an alternative estimator in which we sum over wavelets rather 
than multipoles. 



III. REVIEW OF WAVELET-BASED ESTIMATION 



A. Definition of wavelets 



Wavelets are particularly useful for CMB analysis due to their localization in scale and position. In Rcf. [27; a 
continuous, isotropic wavelet family tp{x,n,R) on was constructed from a 'mother wavelet' ^'(x) by means of 
translations and contractions: ilj{x,n,R) ~ ^'(Ix — n\/R)/R. The mother wavelet should have zero mean and decay 
sufhciently fast at infinity, 



dx ^{x) =0 and 



dx , , < oo . 



(24) 



Each integral is taken over M^. We choose to normalize ^I" so that J dx ^(x)^/i?'^ = 1. 

The wavelet transform of a function /(x) with respect to location n and scale R is defined by w{n,R) — 
J /(x)-0(x, n, R) dx. For a sphere, the location n = n is a unit vector defined by its polar and azimuthal an- 
gles. In this paper we will exclusively use the spherical Mexican-hat wavelet ('SMHW'), 



27r7V(i?) 



1 + 


(1)1 


2 


2 - 


(1)1 













(25) 



where N{R) = R{1 + R'^/2 + R'^/Af/'^ and y = 2tan(6'/2). The SMHW depends only on the polar angle, 9, and the 
scale, R. The Legendre transform of this wavelet, wi{R), satisfies ipsi^^R) = '}2i'Wi{R)Pi{cos 9) [28] . 



B. Wavelets and CMB estimation 



The wavelet transform of a masked CMB temperature map with respect to a set of A'scai scales Ri is 

W(i?„n) = E 

aimWi{Ri)Yi„i{h) , (26) 

hn 

where wi{R) is the Legendre transform of the SMHW. In what follows we redefine the wavelet transform by subtracting 
its mean: we set W{Ri,n) = W(i?i, n) — {\N(Ri, n)), where the average is taken over fi in the unmasked region. This 
has the effect of decorrelating the W{Ri, n) on distances above the resolution scale Ri. It is this property that removes 
the necessity to subtract a linear correction; for further discussion see Refs. ^71 ITS] . 
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Wavelet statistics. The cubic wavelet statistic is defined by 

= ^ /dn WiR,,h)WiRj,n)WiRk,h) , (27) 

where erf = {W{Ri,h)'^), the average again being taken over n in the unmasked region. For isotropic noise we have 
af = (47r)^^ + ^)Ciwf{Ri). In the case of full sky coverage, the expectation value of Wijk is given by 

V^jk = (Wijk) = ^h{Ri)wh{Rj)wi^{Rk)hf^i^ij3i^i^i^ , (28) 



where hi-^^i^i^ is the reduced bispectrum (13) corresponding to the primordial shape whose amplitude we wish to 
constrain; see Eqs. ^ and (|8]). It would ordinarily be computed from the bispectrum predicted by a microscopic 



inflationary model as described in Section IV The quantity hf^^^i^ satisfies 



^2 - (2?i + l)(2^2 + l)(2;3 + l) f\, p , ,p , ,p , ^ ,OQ^ 

hiM3 = ^ J ^ WPh WPh (P) ■ (29) 

Therefore, given the wavelet transform W{Ri,h), the computation of Vijk for a generic bispectrum will involve 



C'(-^scai^niax) Operations. For a real experiment some of the sky must be masked and (28) no longer applies. In this 



case the expectation value of the cubic statistic for each scale must be found using simulations. We describe how 



suitable simulations incorporating the underlying bispectrum can be performed in Section IV 



Optimal estimator. The optimal estimator for the amplitude of the bispectrum shape of interest, bij^i^, i^ [12] 

2-^ijkrst '^3^ ijk,rst ^^st 

where Cijk,rst is the covariance matrix of the cubic statistics Wijk- In terms of the estimator £ introduced in Section [ll| 
it has the schematic form /nl = £/{£)■ We write /^l to indicate that this amplitude depends on the bispectrum 
shape used to construct the estimator. It does not coincide with the traditional definition of Spergel & Komatsu pU] 
unless the bispectrum bi^i^i^ corresponds to the local model. 



In ( 30 ) it is not necessary to include the linear term 

W^r" = 1 / dn W{R,,h){WiR,,h)W{Rk,h)) + 2 perms , (31) 

because of the scale- by-scale subtraction of the mean for each wavelet coefficient P2]- Nevertheless, for completeness. 



we will also subtract this term unless otherwise stated and understand Wijk Wijk — W^jfe"^'^- Eq- (30) can be 
written more succinctly as 



where I,J,K,L are understood as multi-indices, each ranging over a three-component tuple. Therefore / = {i,j,k), 
and similarly for J, K, L. The Fisher estimate for the variance of /nl is 

CTf(/nl) = ./^-i.. • (33) 

To keep numerical errors in the inverse covariance matrix under control we carry out the calculation using principal 
component analysis. Due to the vastly reduced dimensionality — we use of order 10^ cubic statistics — the computation 
is much faster than that of the pixel-by-pixel inverse covariance matrix needed for optimization of bispectrum-based 
estimators. 



In writing this formula we have assumed that the bispectrum bi^i^i^ used to compute Vijk is normalized so that = 1. 



7 



Normalized amplitude estimate 

normalized to the local shape. This redefined measure was introduced in Ref. 
estimator is 



Instead of Eq. (30) we can consider an alternative measure of the amplitude, 

[9|, where it was denoted -Fnl- The 



llijViCjjWj 



1/2 ' 



(34) 



where the Vj°^ should be computed using the local bispcctrum shape. The estimators (30 1 and (34) contain identical 
information and differ only in their normalization. Eq. ( 34 1 is less useful for comparison to specific models because it 
does not yield a correctly- normalized estimate of each amplitude appearing in the primordial bispectrum ([4|. 



IV. PROJECTION FROM PRIMORDIAL TO CMB BISPECTRA AND MAP-MAKING 



A. From primordial to CMB bispectra 

The CMB bispectrum bi-^i^i^ [sometimes called the 'late-time' bispectrum to distinguish it from the primordial 
bispectrum of Eq. ([4|] is a function of the multipoles li, I2 and and therefore may also be decomposed in terms of 
the partial-wave basis Qn- To do so, we define a weighted copy of the bispectrum s(/i, I2, 13) and introduced 'barred' 
coefficients such that 

_ (2^1 + l)V6(2/, + 1)1/6(2/3 -H)i/6 _ 
Shhh = /„ „ „ Khh = a'^Qn[hMM) ■ (35) 

The choice of weighting was explained in Ref. ^ . The barred ('late-time') coefficients should be carefully distinguished 
from the unbarred ('primordial') coefficients which appear in Eq. ( 14 1. For notational simplicity it is sometimes helpful 

(n) 

to renormalize the partial- wave basis by introducing new functions bij^i^ which satisfy 



An) _ ^/CliCl^Cj^ (1 1 1 \ tQR\ 

- + 1)1/6(2^2 + 1)1/6(2/3 + 1)1/6 ^n{il.l2,iz). W 

In terms of this basis we have h^i^i^ = X]„ ^n^hj^is' 

Primordial to late-time mapping. We now aim to express the late-time coefficients in terms of their primordial 
counterparts. To do so we can project sji/aia on to successive basis functions (5n('i, ^2, 'a)- However, the inner 
product ( 16 ) is no longer appropriate because the multipole labels /j are discrete. Therefore we introduce a new inner 
product (see Ref. |5]) defined by 

iJ,g)) ^ KhMM)9{hMM) (2;^ ^ ^^1/3(2^^ + 1)1/3(2/3 + 1)1/3 

'''''' ^ (37) 
j dAi/(Zi,;2,?3)5('i,^2,^3)(2/i + 1)2/^(2/2 + l)'/'(2?3 + l)'/'Ph(M)^b(/^)^i3(M)- 

We first use (20) to express bi^i^i^, and hence si^^i^i^, in terms of the primordial coefficients a^. Projecting the 
resulting expression onto Qm gives 

{{shhi3:Qm{h,l2,l3)} =Ya';^ j Ax X^^nniix) , (38) 

where 

7„™(x) ^ ^ /' Nlll^^{^,,x)Ni-^^{^^,x)N^-JJ^^,x) . (39) 
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In this expression the multi-index n is the triple (ni, 712, 713); the multi-index m is the triple (jni, m2, W3); the bracketed 
indices (■ ■ • ) denote simultaneous symmetrization over both triples with weight unity; and we have defined 



(40) 



Alternatively, projecting the late-time decomposition on the right-hand side of Eq. (35 1 gives 

1 - 

where the same conventions apply for bracketed indices, and 

NnlmA^J■)=^i'2l + ^)^^^Qnl(J^)qrnl(J-'—)Pli^J.), and 7„„ = -!- / d/Z 7V(„^^j (/i)7V„2™2 (/i) JV„3„3) (^) . 



Z 



(42) 



Equating Eqs. (38) and (41) gives the required relationship between primordial and late-time coefficients 

" XI ( / (a;)) 7mn = XI "?rp„ , (43) 



and is defined by this expression. It can be regarded as a projection of the transfer function for the bispectrum 
into mode-space, and describes the change in shape from the primordial era (given by the coefficients a^) to the 
surface of last-scattering (given by the a^). 

The Qn were not specifically constructed to give a good representation of the angular bispectrum, and therefore 
one might harbour some reservations that the approximation of biii2i3 by the same number of basis functions used 
to represent the primordial bispectrum may introduce an unwanted error. However, in practice Eq. (43) proves to 
be extremely accurate, typically producing better than 99% correlation with 0(100) modes. For further details and 
discussion, see Ref. f55l. 



B. Simulating non-Gaussian maps 



It was explained in Section [III| that, for a real experiment, the effects of sky masking and anisotropic noise mean 
that the expectation values Vj — {Wj) for each cubic wavelet statistic must be obtained by numerical simulation. 
[See discussion under Eq. (29).] In Ref. [9;, a simple prescription was given to carry out such simulations. We set 
aim — ^hji + -^NLOtoi' where a^^ is the Gaussian part of each CMB multipole and is a non-Gaussian correction. 



l2m2 hms 

It follows that 



Irn — / , / , ^^IqXz^ vara2ra3 >-» ' V^^y 



a?* a?* 



-,B ^Qn^'-"^ - V V V C^^'a "'2m2"i3m3 

Hm — 2^^n"'lm " g "ll2h^rnm2m3 q q • K^^J 

l2m2 hms ^ ^ 



The afJl^^ can be computed very efficiently, since 

1 Vci, 



B(n) _ -L y'^l 

" 6(22 + 1)1/6 



dn yz„.(n)g(p(V2max)A'C(n)Af,^(n) , (46) 
where the multi-index n is the triple (p, r, s) and the weighted maps Mp{n) are defined as 
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V. APPLICATION OF THE MODAL APPROACH TO WAVELETS 

We are now in a position to connect the modal and wavelet approaches. In particular, we wish to use the modal 
decomposition of some specific primordial bispectrum shape S'^"'^^ (specified by its coefficients a^) to determine the 
expectation value Vi for each cubic wavel et s tatistic. Once these expectation values are determined, the formalism 
of wavelet estimators described in Section 



III 



can be used to recover the amplitude with which slp"^"^ appears in the 

data with a near-optimal error bar. 

To do so, we write Vi = J^n ^n^ni- The matrix Vni can be thought of as a change of basis from partial- waves to 



wavelets and must be computed using the prescription given in Section IV B for evaluation of a non-Gaussian map. 
We fincH 

V^i = — ^ / dn (W^'^(i?,)M^^(i?,)T4^^(")(i?fe) + W''{R,)W'''-^\R,)W''{Rk) + W'''^'^\R,)W'' {R,)W'' {Rk)) , 

(48) 

where the wavelet maps W'^ and are given by (26), with replaced by af^ and afj^"'\ respectively. The 



scale-by-scale mean should be subtracted out as usual. Note that the Vni are independent of any model-specific details, 
such as the shape of the bispectrum for which we are trying to construct an estimator, and can be precomputed. 

Wavelet estimator. It is now possible to write down the wavelet estimator for the amplitude of a bispectrum shape 
specified by primordial coefficients and late-time coefficients a^. It is 

JNL ~ ^ _Q_Q ^-Ijr ■ ^ ' 

l^nm "™ l^IJ ^nlCjj V,nJ 

If desired, this can be written in a form similar to Ref. [5] 

fb ^2n Pn I 

^NL - -Q-Q- ' V^"^ 

where we have defined = VniCyJWj and 7„„ = VniCj^JVmj- 

For any particular experiment, the quantities f3^ and jmn arc model-independent and can be precomputed. (We 
emphasize that they vary between experiments due to the details of noise and masking.) Once these coefficients are 
available, the estimator /^^ fo^' ^''W primordial model can be obtained by trivial summations. 

Orthogonalized modes. Although (|49| is our final result for the wavelet estimator, it can be rewritten in an 
equivalent form which orthogonalizes the partial-wave basis. We perform a Cholesky decomposition of the matrix 
%m to obtain %rn = J2r Kr^ml- Defining = J2n KrO^n and = J2r ^rnP? , it follows that 

This expression is particularly useful because it allows us to deduce that the expectation values (/3^) obtained from 
an ensemble of maps with bispectrum bi-^i^i^ satisfy the relation 

Wn)^o.t (52) 

This relation can be used as a 'sanity check' for simulations of a specific model. 

Alternatively, the analysis could be carried out entirely in wavelet space. Introducing a Cholesky decomposition 
of the inverse covariance matrix CjJ = ^ikLjk we may similarly define = ^kBk / where 



Ak — J2i ^ikVi and — J2j ^JkWj. The 'sanity-check' given by Eq. (52) now becomes {Bx) — Ax 



VI. WMAP7 IMPLEMENTATION 

In Sections [Tl] - [Vl we have assembled the theoretical framework needed to construct wavelet estimators for any chosen 
primordial bispectrum. In this Section and the next we apply this formalism to the coadded V + W foreground-cleaned 



To simplify notation we have omitted the dependence of each wavelet map on n. 
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maps from the WMAP 7-year data release, working up to /max = 1000. The data is at a resolution of 6.9 arcmin, 
corresponding to Nside = 512 for HEALPix [30] . 

Previous analyses of this dataset have used /max = 1500. However, the purpose of this paper is to provide a proof- 
of-concept for the combined modal/wavelet methodology, rather than to obtain the most stringent possible error bar. 
In any case, the bispectrum analysis in Ref. llj, which employed modal techniques, was carried only to /max = 500. 
The authors of that paper noted that the pseudo-optimal approach of Ref. (TU] tends to saturate for larger I. 

For completeness, in Table [l] we list the cosmological parameters used in this analysis. The primordial power 
spectrum is parametrized as a power-law with P<^(k) = A,f,k~'^{k/ki,)"'~^ , and the pivot scale k-i, is taken to be 
k^ = 0.002/1 Mpc"^ 



parameter 


value 




0.0227 


f2c/i'^ 


0.1116 




0.729 


r 


0.085 




1.736 X 10 


ris 


0.963 



Table I: Cosmological parameters used in the WMAP7 analysis. 



A. Generating the wavelets and corresponding masks 



The spherical Mexican-hat wavelet is defined in Fourier space by the Legendre transform of Eq. (25 1 . For each map, 
the wavelet coefficients are obtained by the convolution in Eq. ( 26 ) . We use the same fifteen scales chosen by Curto 
et al. [18], and extend the WMAP KQ75 mask appropriately for each scale. We list the angular scales in Table [Tl| 
together with the fraction of sky available at each scale after applying the mask. 

The construction of an appropriate mask can be done in various ways. Here, our results correspond to masks 
constructed by taking the KQ75 mask without point sources and extending it so that, for the wavelet of scale R, 
any pixel within 2.5i? of a masked pixel is excluded from the analysis. For small-scale wavelets (up to Rq) we 
superpose the mask around point sources. On small scales this mask is believed to be sufficiently extended not to 
cause contamination in the wavelet coefficients. On large scales the effect is negligible. Therefore, further extension 
of the mask for any Ri would be too conservative. We plot the extended masks in Figure [T] 

Constructing extended masks by convolving the existing mask with the wavelet at each scale and leaving out regions 
where the coefficients are contaminated by more than 1% results in very similar masks. 



Wavelet scale 


-Ro 


Ri 


R2 


R3 


R4 


R5 


-Re 


i?7 


Rs 


Rg 


Rio 


Rii 


R12 


Rl3 


Rl4 


Angular scale 





2.9' 


4.5' 


6.9' 


10.6' 


16.3' 


24.9' 


38.3' 


58.7' 


90.1' 


138.3' 


212.3' 


325.8' 


500' 


767.3' 


Sky coverage (%) 


70.6 


70.6 


70.6 


70.6 


70.6 


70.5 


70.4 


70.1 


69.3 


67.3 


63.5 


57.3 


48.4 


36.2 


20.6 



Table II: Proportion of sky covered at each wavelet scale _Ro to -R14. The mask at scale Rq corresponds to the KQ75 mask. 



B. Calculating the modal coefficients 



As we have explained in Section JIB 
coefficients (a^) for the shape 5'^'°'^^ defined in (fTl]) 



r'^ + s'^ and retain only the 80 modes with lowest 



to compute the late-time coefficients (a^) we first compute the primordial 

For each triple (p, r, s) which defines a basis function we compute 
This is sufficient to obtain a correlation of > 99% for 



the usual bispectrum templates (local, constant, equilateral, orthogonal and fiattened), except for the flattened model 
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Ro (unconvolved) 



Ri = 2.9' 



R2 = 4.5' 




Figure 1: (Extended) masks appropriate for the analysis at each wavelet scale. The proportion of sky coverage for each mask is summarized in 



Table 



El 



where we achieve ^ 95% correlation. This can be attributed to an inherent lack of smoothness of the template in the 
flattened limit, which makes the modal decomposition converge only slowly. We use the same number of modes to 
decompose the reduced angular bispectrum hi^i^i^. Later, we will verify that this is sufficient to ensure an accurate 
representation. 

In order to calculate the transfer matrix Tpn from primordial to late-time coefficients we first extract the transfer 

function from CAMB [5S] and compute the line-of-sight projections q^''\x) and q^^'^^{x) defined in Eq. (21). We 
then compute ^mn{x) and 7m„ using Eqs. (40) and (42 1. Combining all these elements enables us to compute the 
transfer matrix from ( 43 ) . 



C. Calculating the observed wavelets and modes 



Our first task is to estimate the expectation value V/ of each cubic wavelet statistic. To do so we simulate 
Gaussian spherical harmonic transforms with the variance for each multipole given by the angular temperature power 
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spectrum, Ci . These are denoted . Using the prescription outhned in Section IV B we generate the simulated non- 



Gaussian multipoles, corresponding to the bispectrum basis function h^j^i^i^ [defined in (36)]. The combination 



aim = off^ + af^''' gives the simulated temperature map. The WMAP7 beam, hi, and noise can be incorporated 
via the transformation 



B{r. 



air. 



aim — hairn + nir 



(53) 



Using Eq. (26) we create simulated Gaussian and non-Gaussian wavelet maps for each scale and apply the appropriate 



masks. Then, for each map, the average in the unmasked region is subtracted. Finally we extract the cubic statistics 
using ( [27| . (Note here that the normalization coefficients for each cubic statistic, trf , are obtained using 
the assumption of isotropic noise. One should not be concerned, because these factors merely represent a normalization 
convention and cancel out in the estimator, Eq. (32).) With the fifteen wavelet scales used in this paper we obtain 680 



cubic statistics. Expectation values for each of the 680 are computed by averaging over 300 simulations, after which 
we compute the change-of-basis matrix Vni using (48). We then compute a wavelet map of the real 7- year WMAP 



data, mask it, and extract cubic statistics Wj in the same way. 

In order to compute the 680 x 680 covariance matrix C/j, we evaluate the expectation value C/j = {Wj^Wf) 



{Wj^){Wj) over 3 x 10"' simulations. The wavelet estimator (32) requires the inverse matrix CjJ, and once this has 
been obtained the quantities and 'Jnm can be computed. 

For each model under consideration we may obtain an estimate for the /nl parameter, f^^, and its expected 



using (32) and (33 1, respectively. We also compute the variance of /^^ from a suite of 



(Fisher) variance, crjp{f^^) 

100 simulations. Irrespective of whether the linear term (31) is subtracted, we recover the expected Fisher variance 
to high accuracy]^ 



D. Validation procedure 



Gaussian validation. To ensure that our implementation is unbiased, we perform 4000 Gaussian simulations. A 
critical 'sanity' check, described by Eq. (52), is to verify that the mean value for each n is consistent with zero, 
within the standard error of the mean. (This is equal to the standard deviation divided by the square root of the 
number of simulations.) We plot the results in Fig. [2] from which we conclude that our methodology successfully 
passes this test. We have evaluated the wavelet-based estimator (51) for these simulations with the result 



nl; 



-0.3 ±23.6. 



(54) 



Note that the standard deviation 23.6 recovers the Fisher value obtained from (33), crF(/NL) = 23.6, precisely. We 
have verified that neglecting the linear term (31 ) results in a minimal (< 2%) difference in the standard deviation. 



Non-Gaussian validation. We produce 100 local simulations with /jlj^ = 100, using the method described in |31j . 
For each of these simulations we extract the observed modes, /3^. The critical sanity; check (52) can now be carried 



out by comparing the theoretical expectation a„ to each observed mode. In Figure 
satisfied within the error bars. The recovered value of /nl is found to be 

(/nl) = 99.9 ± 2.5 , 
where the error bar quoted in this case is the standard error of the mean. 



we show that this test is again 



(55) 



VII. 7- YEAR WMAP CONSTRAINTS 



In this Section we use the methodology described in Section [Vl| to obtain constraints on a selection of nearly scale- 
invariant models. To date, most bispectrum analyses have considered the local and equilateral models, owing to their 
physical significance and computational simplicity. We extend these to include the constant, orthogonal and flattened 
templates. In Ref. ^llj these models were studied using modal methods and a bispectrum-based analysis. However, 



We replicate the result of [17], finding that the variance without accounting for the linear term is within < 2% of the variance when this 
contribution is included. 
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0.06 



0.02 




-0.02 



-0.06 



Figure 2: Mean of (/3^ ) from 1000 Gaussian simulations. To verify eonsistency of eaeh mode with zero we plot the standard error of the mean 
(ie. the standard error divided by the square root of the number of simulations). Every mode is eonsistent with zero within two stan- 
dard errors of the mean. 




40 

Mode Number 



Figure 3: Comparison of the modal eoefficients for the loeal shape, , to the observed coefficients, . The observed coefficients are com- 
puted by simulating 100 local maps w ith — 100 as described in Ref. [31 . We plot the mean of these modes and the standard 
error of the mean. The 'sanity' check ||52[ is satisfied. 
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the wavelet-based approach adopted here yields a treatment which is closer to optimal. For example, in the case of 



local non-Gaussianity the error bar is reduced from A/nl = 27.6 to A/nl = 23.6. In Table III we compare our results 
with those of |llj . highlighting the improvement in optimality achievable with the approach adopted in this paper. 
The comparison highlights the significant improvement over a standard modal-based approach, which, thus far, has 
been constrained in its scope, to using the same basis for the estimator as that used for the decomposition of the 



shape. This necessity may be avoided by employing the change of basis matrix, ( 48 ) , with the modal technique used 



for the shape decomposition and wavelets being employed for extraction of the data in this work. 



Shape 


Current Paper 


FLS XL 


Local 


38.4 ±23.6 


20.3 ±27.6 


Constant 


-10.1 ±60.6 


30.5 ±95.9 


Equilateral 


-119.2 ± 123.6 


1.9 ± 127.4 


DBI 


-50.1 ± 104.2 


17.1 ± 121.8 


Orthogonal 


-173.2 ± 101.5 


-51.4 ± 103.8 


Flattened 


6.6 ± 10.4 


3.0 ± 10.9 



Table III: Comparison of the WMAP7 constraints found in this worlt with those of Fergusson, Liguori and Shellard 1111 . The improvement is 

acliiovcd by employing a more optimal estimation method, allowed due to the freedom in making a different choice of basis to extract 
the data from that used to decomposed the theoretical shape. 



A. Local model 

The local model is defined by a series expansion of the primordial gravitation potential, $, in powers of an exactly 
Gaussian potential, <f>G- It gives [35] 

$(x) = $g(x) ± /nl($g(x)2 - ($g(x)2)) . (56) 
Neglecting scale dependence, the local model is an accurate match for the non-Gaussian contribution generated by 



interactions which operate on superhorizon scales. It generates a bispectrum given by ( 10 ). Examples of scenarios that 



may produce appreciable local non-Gaussianity include multifield inflation and curvaton scenarios. By construction 
the shape S^°'^^ — 1 for this model. For a comprehensive review, see Chen [5] and references therein. 

In Figure [5] we plot the primordial bispectrum shape. The dominant signal occurs in the corners of the triangle, 
corresponding to 'squeezed' configurations where one momentum is much smaller than the other two. When transferred 
to the CMB bispectrum, this (almost) scale-invariant shape is redistributed, resulting in the presence of peaks in the 
three-dimensional CMB bispectrum. However, the dominant signal remains along the edges of the tetrahedral domain, 
i.e. where one I is much smaller than the other two. 

In Figure |4j we compare the modal coefficients for the local model against the coefficients reconstructed from 7- 
year WMAP data. We also plot the cumulative sum J2n=cr '^nPn / Yl'n=Q^'^nY to establish that the estimator does 
converge with 80 modes. Our final constraint on the amplitude of a local-type bispectrum is 

/nl = 38.4 ± 23.6 or F]^1 = 38.4 ± 23.6. (57) 

This result is competitive with the outcome of other wavelet-based analyses. For example, Donzelli et al. |T7] quoted 
the constraint /^l = 37. 5 ±22.3. The small difference in our results may be explained by the use of a slightly different 
masking procedure, and (perhaps more importantly) because their analysis used data up to ^max = 1500. 

B. Constant model 

The constant model is simply S'$(fci, fc2, /ca) = 1. It is interesting because it produces a CMB bispectrum due 
completely to the transfer functions. One possible microphysical realization may occur during an epoch of quasi- 
single field inflation [33J. Our wavelet-based estimator yields the constraint 

jco^nst ^ ^ _^ QQ g ^co^nst ^ _3 g _^ 23.6. (58) 
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_6 I 1 1 1 1 1 1 1 1 

10 20 30 40 50 60 70 80 

Mode Number 

(a) Comparison of the modal coefficients of the local model against the recovered modes of the 
WMAP bispcctrum. The coefficients are normalised to aid comparison. 




.20 I 1 1 1 1 1 1 1 1 

10 20 30 40 50 60 70 80 

Maximum Mode Number 

(b) Cumulative contribution of each mode to the amplitude /nl- We plot the partial sum 

^n=o^ X]n=oC^n )^ agalust the maximum mode number A^max- Despite the slow con- 

vergence of the WMAP signal, the estimator for /nl has converged. 

Figure 4: Modal comparison for local model. 



C. Equilateral and DBI models 

Interactions which operate on supcrhorizon scales produce a local-shape bispectrum because causality requires the 
interaction to consist of a long-wavelength modulation of the background experienced by the short modes. This 
correlation between long and short modes is maximized in the squeezed limit. 
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Figure 5: Canonical shape function \12\ for the primordial local bispeetrum. Due to the (almost) scale invariance of the shape it is only nec- 
essary to plot a particular slice. The parametrization {a,P) for each slice is chosen as described in Ref. [9]. In particular, k\/kt — 
(1 + a + fca/fct = (1 - a + 0)/4, fca/fct = (1 - /3)/2, whore kt = Ei k^. 



In comparison, interactions which dominate on subhorizon scales typicaUy produce no signal in the bispeetrum, 
because subhorizon modes fluctuate incoherently and average to zero. An exception, where the initial state is non- 
empty, will be considered in Section p/II E| below. Neglecting that possibility, significant effects can be produced only 
near the epoch of horizon exit, where the fluctuations are beginning to behave coherently. In canonical slow-roll, 
single-field inflation the interference between horizon-scale fluctuations does not generate significant non-Gaussianity. 
However, with non-standard kinetic terms the amplitude of these effects may be enhanced [34-36 . Examples include 
DBI infiation and /c-inflation. 

The equilateral template is a separable approximation to the bispeetrum produced by such models. It produces 
strong correlations for roughly equal k because it is dominated by interference effects between wavenumbers which 
leave the horizon nearly simultaneously. We plot the shape function for the DBI model and the equilateral template 
in Fig. [6j They correspond to 



rDBI 



3 7„2n 



P<i.(fci)P<i,(fc2) + 2 perms - 2 



{kfkjki - Akfk]ki) 

2 /3 

P*(A:i)F$(A:2)F$(fc3)' 



(59) 



Pl'\ki)P;,'%k2)P^{k:i) + 5 perms 



,2/3 



(60) 



Using (19) it can be shown that correlation between these shapes is 98%. The wavelet-based estimator gives the 



constraints 



j-DBI _ 
/NL - 

•/NL 



-50.1 ±104.2 
-119.2 ± 123.6 



pDBI 



-11.4 ±23.6, 
-22.8 ±23.6. 



(61) 
(62) 



A variety of other bispectra dominated by interference effects near horizon-crossing, including the case of ghost 
inflation, were considered in Ref. [IT. All these models are highly correlated with the equilateral template and 
produce similar constraints. 



D. Orthogonal Model 



irth 



The orthogonal shape is a linear combination of the constant and equilateral shapes, corresponding to S[ 

(2/3)5*1°"''* lUISZ]- It is roughly orthogonal to both the equilateral and local shapes. A separable template for 
its bispeetrum is 



Qcq 



b: 



lOrth 



6 3 



Pl'\k^)Pl'\k2)P^{k:,) + 5 perms 



P^{ki)Pi{k2)Pi{k3 



2/3 



loc 



(63) 
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DBI shape S^^ 




Equilateral shape 



equil 




(a) DBI shape 



(b) Equilateral shape 



Figure 6: Canonical shape function, jl2^ , for the primordial DBI bispcctrum and equilateral bispeetrum. The high degree of correlation be- 
tween the two shapes, apparent from these plots, is 98%. 

It has been shown to arise in the DBI-Gahleon model by Renaux-Petel ; see also Ref . |,39 . We find the following 
constraint 



Ji 



■orth 
NL 



-173.2 ±101.5 or F^'^^ 



-40.3 ±23.6. 



(64) 



We note the consistency of this result with the constraint /nl'^ = —159.8 ± 115.1 obtained by Curto et al. [T5] . 



E. Flattened Model 



If the initial state for fluctuations is not empty then it is possible to produce a bispeetrum describing maximum 
correlation for 'flattened' configurations where, eg., ki sa A:2 + ^3- The correlation arises because, with a nontrivial 
initial state, it is possible to find a 'negative' energy fluctuation with time dependence ~ e"''^^* which interacts 
coherently with two positive energy fluctuations with time dependence ~ g"'"'''^*, e+''^^*. When fci w ^2 + fca the 
interaction is coherent over arbitrarily long times and does not average to zero in the subhorizon era. 

Holman & ToUey studied a model in which the bispeetrum produced by this effect was [JU] 



B. 



flat 



6 



39(fciA:2fc3)2 



"-2 "-3 



^2^3 



2 perms +12 + 8 



kik2 + kikz - k2kz 
(/c2 + fc3-fci)2 



2 perms 



(65) 



This is not separable. Its analysis is computationally intensive without a method such as modal decomposition, 
although alternatives approaches exist such as the use of Schwinger parameters |41j . 

Eq. (65) diverges in the flattened limit because the interaction continues over arbitrarily long times, and therefore 
becomes sensitive to whatever physics was operative throughout the inflationary era. By comparison, inflationary 
predictions using a vacuum initial state decouple from this unknown high-energy physics. To handle the divergence we 
parametrize our ignorance of the relevant physics using a cutoff [9], setting the bispeetrum to zero for fci + fc2 — fca < Zkt 
(or its permutations), where the perimeter kt = fci + fc2 + fc3 was defined below Eq. ( 15 ). In this paper we take Z — 0.03 
as a fiducial value, although in a dedicated analysis Z should be allowed to float. Furthermore, we smoothen the 
shape near the edges by employing a low pass (Gaussian) filter. With this choice, our constraints on the flattened 
model from the 7-year WMAP data are 



I- flat 

/nl 



6.6 ± 10.4 or F. 



fiat 
NL 



15 ±23.6. 



(66) 



VIII. CONCLUSIONS 



In this paper we have developed a framework which combines partial-wave or 'modal' techniques with wavelet-based 
estimators for the CMB bispeetrum. 
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Wavelet-based techniques are particularly efficient for CMB analysis because they take advantage of simultaneous 
localization on the temperature map in both scale and position. However, it is not straightforward to build a 
wavelet-based estimator for the amplitude of an arbitrary primordial bispectrum _B$. Combining the wavelet-based 
methodology with the decomposition of i3$ into a basis of partial-waves enables an efficient analysis — whether i?$ is 
separable or not, provided only that it is a relatively smooth function of wavenumber. 

Our framework has other advantages. Optimal bispectrum-based estimators are typically hampered by the need 
to invert a pixel- by-pixel covariance matrix of size ^ 10^ x 10^. Because the wavelet-by- wavelet covariance matrix is 
typically much smaller, of order ~ 10^ x 10"^, the problem is numerically much more tractable. Despite this large gain 
in numerical efficiency there is comparatively little trade-off in optimality. Indeed, our final error bars are quite close 
to those achieved by the optimal pixel-by-pixel approach. In addition, due to their localization properties, wavelets 
have been shown to allow for accurate analysis of point-sources, foregrounds, and other systematics [111 132]- Hence, 
integrating the partial-wave methodology with a wavelet-based analysis has the potential to reproduce the successes 
of both. The work presented in this paper generalizes the scope of wavelet-based estimation to allow for analysis of 
arbitrary primordial bispectra. 

We have implemented our methodology for the 7-year WMAP data. Our constraints are competitive (to within 
~ 5 — 10%) with comparable constraints published elsewhere, and represent an improvement of up to ~ 15% in 
comparison with the bispectrum-based modal estimator of Ref. [TT]. In any case there is much to be gained from 
implementing different estimators: they will be sensitive to different combinations of the data, including the underlying 
systematics. 

In future work, we intend to study the efficacy of the method in more detail and pursue an extension to the 
trispectrum. (See also Refs. [1TJ[52].) Although our constraints show that the 7- year WMAP data are consistent 
with Gaussianity, we will shortly be presented with an improved data set from Planck. We hope that the techniques 
described in this paper can help to correctly categorize the source of any signal — whether of primordial origin, or a 
foreground. 
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